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Abstract 

The cross-over from coherent to incoherent exciton transport in disordered polymer films is 
studied by computationally solving a modified form of the Redfield equation for the exciton density 
matrix. This theory models quantum mechanical (ballistic) and incoherent (diffusive) transport 
as limiting cases. It also reproduces Forster transport for certain parameter regimes. Using model 
parameters appropriate to polymer thin films it is shown that short-time quantum mechanical 
coherence increases the exciton diffusion length. It also causes rapid initial energy relaxation and 
larger line widths. The route to equilibrium is, however, more questionable, as the equilibrium 
populations of the model do not satisfy the Boltzmann distributions over the site energies. The 
Redfield equation for the dimer is solved exactly to provide insight into the numerical results. 

PACS numbers: 71.35.-y, 78.20.Bh, 78.55.Kz 
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I. INTRODUCTION 



Energy transport via exciton migration in disordered organic semiconducting solids has 
been extensively investigated by theoretical and computational modellin g 1 ! 2 ' 3 ! 4 ' 5 . It is gen- 
erally assumed that the dephasing rates in these systems are large compared to the exciton 
transfer rates, so that the quantum mechanical phase memory of the exciton is rapidly 
lost, rendering the exciton motion incoherent. The effective medium approximation was 
developed by Movaghar and co-workers^ in the incoherent limit to investigate the role of 
spatial and energetic disorder on charge and energy transport properties. This mean-field 
theory agrees remarkably well with Monte Carlo simulations at high enough temperatures, 
predicting a roughly — log(i) time dependence for the energy relaxation. At low tempera- 
tures, however, mean-field theory overestimates the relaxation pathways 4 -, and thus fails to 
predict the 'freezing-in' of the energy relaxation observed both experimentally and in Monte 
Carlo simulations. For a Gaussian distribution of energetic disorder the pseudo-equilibrium 
diffusion coefficient, D(T), is found to behave as, 



at high temperatures (where T is proportional to the width of the Gaussian distribution), 
while at low temperatures it is exponentially activated— 

Conversely, it is generally assumed that dephasing rates in biological light harvesting 
complexes are generally small enough that the exciton motion remains coherent for times 
long enough for the exciton to successfully reach the reaction center before it can radiatively 
recombine. Indeed, in light harvesting complexes the ballistic nature of exciton transport 
coupled to an energy landscape that 'funnels' the exciton to the reaction center leads to 
particularly efficient energy transport.-^ In a disordered energy landscape with sufficiently 
large energetic disorder, on the other hand, localization occurs in the coherent regime^. 

Dephasing rates in typical organic compounds are significantly enhanced by intra and 
inter molecular disorder that scatters the exciton wavefunction 9 . Thus, typical dephasing 
times are roughly 1 ps in an ordered conjugated polymer chain^, compared to roughly 100 fs 
in a disordered system^. Nonetheless, the advent of ultra-fast spectroscopy allows relaxation 
phenomena occurring on sub- 100 fs time-scales to be observed, implying that theoretical 
modelling of exciton transport should also take into account quantum coherence effects on 
these time-scales. Indeed, more recent Monte Carlo simulations of exciton transport in the 




(1) 



2 



incoherent limit fail to predict the observed very fast relaxation processes that lead to a 
rapid reduction of the average energy and a broadening of the spectral lines 5 . Although this 
rapid energy reduction was attributed to vibrational relaxation 5 -, as will be demonstrated 
in this paper, another possible mechanism is coherent energy transport. The 'anomalous' 
red-shifted emission also emphasizes the role of quantum coherence on exciton dynamics, 
as this has been attributed to the recombination of an exciton delocalized over a pair of 
polymer chains 5 -. 

As for light harvesting complexes, the successful operation of a polymer photovoltaic 
device requires that a photo-excited exciton migrates to a polymer heterojunction and then 
disassociates before it can recombine. An understanding of the role of quantum coherence 
and energetic disorder on exciton transport is therefore essential if realistic simulations of 
device properties are to be performed. 

This paper presents computational studies of the cross-over from coherent to incoherent 
exciton transport. The model that we adopt assumes a localized basis whereby an exciton 
initially created on a chromophore migrates to neighboring chromophores. This model re- 
tains the properties of fully coherent and fully incoherent transport as limiting cases. In 
particular, for short times the quantum mechanical processes lead to a coherent wavepacket 
that is delocalized over more than one chromophore. At long times dephasing processes de- 
stroy phase coherence, and the motion becomes 'classical'. Forster-like dipole-dipole induced 
exciton transport is described by a limiting case of this model. We show that short-time 
coherent motion leads to increased exciton diffusion lengths, and to initial rapid energy 
relaxation and increased line widths. 

The plan of this paper is as follows. The next section outlines the model, while Section 
III describes the numerical techniques to solve it. Section IV describes and discusses our 
results, and we conclude in Section V. 

II. THE MODEL 

Exciton transfer from chromophore m to chromophore n, mediated by the dipole-dipole 
interaction, is parameterized by the energy transfer integral, J' mn . Assuming that each 
chromophore exists either in its ground state or a single excited state, and that the excitation 
energy of the mth chromophore is e m , the Hamiltonian that describes the coherent exciton 
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dynamics in the absence of exciton interactions is, 

H= ^e m \m)(m\ + ^2(J' mn \m)(n\ + J' nm \n)(m\) . (2) 

m inn 

The basis state \m) represents an exciton localized on the mth chromophore. In general, 
exciton interactions, particularly with defects and the heat bath via the exchange of phonons, 
will cause a loss of phase coherence. This behavior is conveniently modelled by an equation 
of motion for the reduced density operator, p, defined by 

p = Tr{W(t)}, (3) 

where W(t) is the full density operator and the trace is over all the degrees of freedom of 
the environment. 

In a localized exciton basis the matrix elements of the reduced density operator are 
Pmn = (m\p\n). An equation of motion for these matrix elements may be formally de- 
rived by performing the trace in Eq. (JH|— . In this paper, however, we assume a semi- 
phenomenological approach and adopt a Redfield-like equation that describes coherent and 
incoherent processes as limiting cases, and in general models the cross-over from short-time 
coherent behavior to long-time incoherent behavior. This model is defined as follows by 
decomposing the equation of motion for the matrix elements into its constituent parts: 

dpmn{t) 



dt 



L\Pmn{t) + L 2 p mn (t) + L 3 p mn (t) + L 4 ft mn (t), (4) 



where 



LlPmn — —i [H, p] mn — —i^mnPmn ~ T {JmlPln ~ JlnPmi) i (5) 



h 



J->2pmn $mn ^ ^ {J^mlPmm - ktmpu) > (6) 

e 

L?>pmn 2^y mn (l $mn) Pmni (7) 

■LiPmn $mnPmn / I'm- (8) 

In these equations we have defined, J mn = J mn /h, u mn = (e m -e n )/h, 2^ mn = 7 m + 7n + 27o, 
and 7 m = |X)m-^*W- 



and 
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L\p represents the coherent, ballistic motion of the exciton, and describes the exciton 
motion in accordance with the time-dependent Schrodinger equation, 

it£^=Hm)), (9) 

where H is given by Eq. (J2J). 

L2P represents the incoherent, diffusive motion of the exciton arising from population 
transfer from chromophore to chromophore. Formally, this is associated with vibrational- 
induced exciton transfer via the spatial modulation of the dipole-dipole coupling^. With 
this term alone Eq. (0J) is equivalent to the Pauli master equation, 

= - Yl ( k ^ p ™ - k ^ p t) » ( lo ) 

t 

where P m = p m m- Energy relaxation occurs when the rate for energy transfer to a higher en- 
ergy chromophore is smaller than the rate for energy transfer to a lower energy chromophore. 
It is customary to assume that the rates satisfy, 

^ = exp(-hLu em /k B T), (11) 

which guarantees that the equilibrium populations satisfy the Boltzmann distribution for 
Eq. (jlOjl . As will be shown later, however, Eq. (flT]) does not guarantee that the equilibrium 
populations of the Redfield equation (Eq. (J1J) satisfy the Boltzmann distribution. 

L 3 p represents the coherent (or transverse) dephasing of the off-diagonal elements of the 
density matrix. In the Bloch model the transverse dephasing time, T 2 = (27) ~ 1 . 

Finally, L^p represents the population decay via exciton recombination mechanisms. If 
the dominant decay is via radiative recombination then, 

Tm = r (—) , (12) 

where e is the average chromophore energy. 

Eq. (j3J) has been extensively studied for translationally invariant systems (se o 1 ?, 14 oi"£, 
for example). Exact results for the mean-square-displacement, (r 2 ), have been obtained by 
Reineker— for a c/-dimensional cubic hyper-lattice. For an exciton created at the origin at 
time t — 0, 

<r 2 (t)> = a 2 d ^(k + jjt + ^ (exp(-27t) - 1)^ , (13) 



where a is the lattice parameter. For long times > 7 x ) this expression reproduces the 
random walk result, 

(r 2 ) = 2a 2 d + t + constant, (14) 

with an additional Forster contribution to the diffusion (namely, the J 2 /7 term) arising 
from the dipole-dipole coupling. This result reflects the fact that Eq. (j3J) approximates to 
the diffusion equation in the long-time limi^. Thus, in this model in the long-time limit 
the dipole-dipole interaction contributes to incoherent motion in two ways. First, explicitly 
from LiPmn (where k formally is proportional to the spatial derivative of J mri — ) and second, 
implicitly from the Forster mechanism. 

For short times (t — > 0), the expected ballistic result is modified by an anomalous linear 
term resulting from the assumptions in the derivation of Eq. ()13j) 13 : 

(r 2 ) = 2a 2 d (kt + J 2 t 2 ) . (15) 

It is immediately apparent from Eq. (fT3|) that for short times (but for t > k/J 2 ) the exciton 
transport is dominated by coherent, ballistic processes, and therefore the exciton travels 
much further than it would if only incoherent, diffusive processes are considered. 

Finally, the fully quantum mechanical behavior of Eq. Q is also reproduced from Eq. 
(Unj) by setting k = 7 = 0, 

(r 2 ) = 2a 2 dJ 2 t 2 . (16) 

Appendix [X] describes the solution of Eq. for the special case of a dimer. For equal 
site energies this solution illustrates the decay of coherence and the Forster limit when 
7 3> J, while for unequal site energies it illustrates the persistence of coherence and that 
the equilibrium populations are not determined by the Boltzmann distribution over the site 
energies when the rates satisfy Eq. (jllj) . To our knowledge no analytical results exist for 
the transport properties associated with Eq. (|IJ) in a general disordered energy landscape. 
The numerical techniques to solve this model are therefore described in the next section. 



III. METHODOLOGY 

Eq. (HJ) may be expressed as the set of M x M simultaneous equations, 



where M = iV s 2 ites and iV S i tes is the number of sites in the lattice, yi(t) = p mn (t), and A is 
the coupling matrix. The formal solution of Eq. (fT7|) is, 

Vi {t) = y £S iJ exp{\ j t)S£y k {0). (18) 

jk 

Here, S is the matrix whose columns are the eigenvectors of A, {A} are the eigenvalues of 
A, and {y(0)} are the initial conditions. Assuming that the mth chromophore is excited at 
t = 0, we have, 

Vi(t) = Y;S id exp(\ j t)S£. (19) 

j 

The evaluation of Eq. (|T9"j) requires a diagonalization of the M x M sparse, complex, and 
non-hermitian transformation matrix, A. For most problems of interest this matrix is far 
too large to numerically diagonalize completely. However, the real parts of the eigenvalues 
of A are the rates for the decay of the eigenmodes. Consequently, for times greater than 
an arbitrary cut-off, r, only eigenvalues whose real parts satisfy 9ft|Aj| < r _1 need to be 
computed. If only this long-time behavior is required, it is only necessary to diagonalize A 
for a small sub-set of the entire spectrum using sparse-matrix diagonalization techniques^. 
This approach is further described in Appendix [Bj The long-time behavior can then be 
matched to the short-time behavior obtained by standard numerical integration techniques, 
as described below. 

In practice, the equation of motion for the density matrix, Eq. (J%J). may be solved more 
efficiently for large systems by standard numerical time-discretization techniques. To ensure 
both stability and accuracy the Crank-Nicolson scheme - an average of implicit and explicit 
forward-time-centered-space discretization methods - is employe d 15 ^ 6 . In order to adapt the 
Crank-Nicolson scheme for Eq. (JIJ it is necessary to adopt the Operator Splitting Method 17 
by decomposing the spatial differential operator on the right-hand-side of Eq. (J3J) into a sum 
of zero or one-dimensional operators. 

Expectation values of operators corresponding to dynamical variables are found from the 
density matrix in the usual manner, via, 

(0(t)) = J2prnn(t)O nm . (20) 

mn 
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IV. RESULTS AND DISCUSSION 



For an iV S ite site lattice the number of matrix elements (or components) of Eq. (jlj) is Nf iiB , 
rendering this a computationally very expensive problem. To obtain numerical results from 
time-scales of 10~ 15 to 10~ 9 seconds it has therefore been necessary to restrict the size and 
dimensionality of the lattice. In this paper we describe results for two-dimensional square 
lattices of up to 41 x 41 sites. 

Our parametrization of the model follows closely that of Meskers et a/.—. We take nearest 
neighbor interactions on a square lattice, with k = 0.00088 eV (= (5 ps) -1 ), r = 300 
ps, 70 = 0.0207 eV (= (200 fs) -1 ), while J (which models the strength of the quantum 
coherence) is an independent parameter, varying from 0.02 eV to 0.00414 eV (= (1 ps) -1 )^. 
The value of J = 0.00414 eV corresponds to a Forster transfer rate of </ 2 /7o = 0.00088 
eV, so in the limit that k = and long times the model maps directly onto the Mesker 
parametrization. The energetic disorder is modelled by a Gaussian distribution function^ 
of mean energy, e = 3.0 eV and a = 0.07 eV. We perform averages over ensembles of five 
realizations of the disorder. The system is initially excited at the origin at t — with an 
energy of 3.1 eV. 

A. Root-mean-square distance 

We first investigate the role of quantum coherence on the root-mean-square distance, 
(r 2 (t)) 1//2 , travelled by the exciton at a time t. Fig. ^ shows (r 2 ) 1 / 2 for various values of J. 
As expected, larger values of J lead to greater distances travelled. For a value of J = 0.02 
eV the exciton has reached the boundaries of the 41 x 41 square lattice within its mean life 
time, tq = 300 ps. In contrast, in the classical limit (J = 0) the exciton has only travelled 
ca. 4 lattice units in this time. 

We next discuss the role of energetic disorder on the root-mean-square distance. As 
Eq. indicates, in the incoherent limit the scaled pseudo-equilibrium constant, D/D , 
is a function of disorder (and temperature) only. In contrast, for purely coherent motion 
diffusion is related to both the exciton band width and the scale of the disorder, and indeed 
localization occurs for sufficiently larger disorder—. Eq. (|13|) indicates that (r 2 ) 1 / 2 / J is a 
function of 7 and k/J 2 = k'. Thus, in the absence of disorder plots of (r 2 ) 1 / 2 /J versus 
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FIG. 1: The mean root-mean-square distance, (r 2 ) 1//2 , versus time for various values of J and 
lattice size: J = 0.02 eV (solid curves); J = 0.01 eV (dotted curves); J = 0.00414 eV (short- 
dashed curves); J = (long-dashed curves); 41 x 41 lattice (no symbols); 25 x 25 lattice (x). Also 
shown is the mean recombination life time, To = 300 ps. 

time for fixed values of 7 and k' will coincide for all values of J. Fig. |2] shows that, as 
expected, disorder reduces the value of (r 2 ) 1//2 . Rather unexpectedly the scaled root-mean- 
square distance increases with decreasing J. This is an artefact, however, of the fact that 
disorder is more effective at hindering transport when quasi-equilibrium is reached, which 
takes longer to achieve for smaller values of J. A full analysis of the scaling of diffusion with 
disorder requires further more extensive numerical calculations and analytical studies. 

B. Energy relaxation 

The expectation value of the energy as a function of time is obtained via Eq. (j20j) with 
O = H, and H given by Eq. (J2J). As shown in Fig. El the initial rate of energy relaxation 
increases as J increases. By 10 -12 s the energy has relaxed by 0.1 eV when J = 0.02 eV and 
by 0.05 eV in the classical limit (J = 0). This result can be understood by the observation 
that in the coherent regime the exciton initially forms a wavepacket that travels ballistically, 
and hence rapidly samples a wide ensemble of sites. 

Notice, however, that the equilibrium value of (E) increases with increasing J, and for 
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FIG. 2: Scaled root-mean-square distance, (r 2 ) 1 / 2 /J, versus time for various values of energetic 
disorder and k' = k/J 2 on a 41 x 41 lattice, k' = 2.2 eV" 1 (solid curves); k' = 51.3 eV _1 (dotted 
curves); a = 0.07 eV (curves with symbols); a = (curves with no symbols). The exact results 
are obtained from Eq. (|13|) . 

non-zero J it obviously does not satisfy the canonical ensemble value given by the classical 
(Boltzmann) distribution over the site energies {e m }. This result is a consequence of the 
fact that for this model with a disordered energy landscape the equilibrium values of the 
coherences (i.e. the off-diagonal density matrix elements) are not identically zero. Thus, 
because of quantum mechanical derealization, the probability for occupying a high energy 
site is higher than that predicted by the Boltzmann distribution. This increase in the poten- 
tial energy caused by the derealization onto higher energy sites is not entirely compensated 
by the kinetic energy reduction arising from the quantum mechanical derealization. A full 
treatment of a dimer with unequal energies is given in Appendix El to further illustrate this 
point. 

Fig.Hshows the line width, (AE) = (E 2 - (E) 2 ) 1 / 2 , for various values of J. As expected 
from the previous discussion concerning the initial faster derealization of the exciton in the 
quantum mechanical limit, the line width increases more rapidly for larger values of J. The 
equilibrium line widths are also larger for larger values of J, because the measurement of 
energy is with respect to the Hamiltonian, Eq. (|2|L whereas the system is not in an eigenstate 
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FIG. 3: The mean energy, (E), versus time for various values of J for the 41 x 41 lattice: J = 0.02 
eV (solid curves); J = 0.01 eV (dotted curves); J = 0.00414 eV (short-dashed curves); J = 
(long-dashed curves). 

of that Hamiltonian. 



V. CONCLUSIONS 



The successful operation of a polymer photovoltaic device requires that a photo-excited 
exciton travels to a polymer hetero junction and then disassociates before it can recombine. 
However, exciton transport in molecular systems is strongly dependent on molecular order, 
as this both determines the exciton dephasing times and the energy landscape through which 
the exciton travels. 

This paper describes numerical solutions of a modified Redfield equation for the exciton 
reduced density matrix in order to study the role of quantum mechanical coherence and 
energetic disorder both on the exciton diffusion length and its energetic relaxation. Using 
model parameters appropriate to polymer thin films, we showed that increased quantum co- 
herence (achieved by increasing the exciton band width) leads to increased exciton diffusion 
lengths. It also leads to initially more rapid energy relaxation and to wider line widths, in 
qualitative agreement with experiment^. Quite generally, increased disorder implies shorter 
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FIG. 4: The mean standard deviation of the energy, {AE), versus time for various values of J 
for the 25 x 25 lattice: J = 0.02 eV (solid curves); J = 0.01 eV (dotted curves); J = 0.00414 eV 
(short-dashed curves); J = (long-dashed curves). 

coherence times and energetic trapping, thereby strongly inhibiting the excitons ability to 
migrate successfully to a polymer heteroj unction. 

Although the model is appropriate for the short-time transport and energy relaxation 
processes, its applicability for describing the long-time route to equilibrium is more ques- 
tionable. In particular, we showed that the particular assumption of Eq. pip for the rates 
for population transfer do not in equilibrium reproduce the Boltzmann distribution over the 
site energies of Eq. (J2J). (This result was proved for the special case of the dimer.) The 
cause of this discrepancy is the choice of basis for the exciton transport. As described in 
re£i£, taking as a basis the exciton eigenstates of Eq. (J2J) would give equilibrium described by 
the Boltzmann distribution over these eigenstates. Unfortunately, however, this prescription 
does not correspond to the physically intuitive picture of local exciton transfer from site to 
site (and consequently does not reproduce Forster transfer as a limiting case). 

Further work will include simulations over larger lattices in three dimensions and the 
possible development of a mean field theory in order to understand the relations between 
disorder and quasi-equilibrium diffusion. 
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APPENDIX A: EXACT SOLUTION OF THE REDFIELD EQUATION FOR THE 
DIMER 



The exact solution of the Redfield equation of motion of the density matrix for the dimer 
illustrates a number of important results. For equal site energies it illustrates the decay 
of coherence and the Forster limit when 7 ^> J. For unequal site energies it illustrates 
the persistence of coherence and that the equilibrium populations are not determined by 
the Boltzmann distribution over the site energies when the rates satisfy the expression, Eq. 

The equation of motion of the matrix elements for the dimer in the absence of population 
loss are, 

11 : -iJ{p2i - P12) - k n pu + hiP22, (Al) 



and 



dt 
dt 

9p22 

dt 
dp 2 i 



-iup 12 - iJ (p 22 - P11) - 270P12, (A2) 
-- iJ (P21 - P12) + fciaPii - ^2iP22, (A3) 
= iup2i + iJ (P22 - P11) - 270P12, (A4) 



dt 

where u> = (e x — e 2 )/h and we have set 7 = j . The initial conditions are taken as pn(0) = 1 
and p 22 (0) = pi 2 (0) =p 2 i(0) = 0. 

When the site energies are equal ku = k 2 i = k, u = 0, and the equations are easily solved 
by Laplace transforms to give, 

P l {t) = p n (t) = ~ (l + (cos(Qt) + (2^*j sin(m) J exp(-(A; + j )t) J , (A5) 



P 2 (t) =P22{t) = 1-Pi, (A6) 

K( Pl2 (i)) = 3?(p2i(i)) = 0, (A7) 

and 

2 7 

5(pi 2 (t)) = -^(P2i(t)) = ■^-sin(nt)exp(-(fc + 7o)*), (A8) 
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where 

0= (4J 2 -(A;- 7 o) 2 ) 1/2 . (A9) 

Notice that the real components of the coherences are identically zero. Also, for under- 
damped systems the populations decay to their classical values of 1/2 and the imaginary 
components of the coherences decay to zero on a time-scale of (k + 7o)~ 1 - 
The Forster limit is derived for the case that k = and 70 3> J. Then, 

Pi(t) - \ (l + exp(-2fct)) + 0(J/ 7 o) 2 , (A10) 

where k, the Forster rate, is k = J 2 /7o, and Pi is a solution of the Pauli master equation, 
Eq. (fTU|) . with k£ m replaced by k. 

When the site energies are unequal the transfer rates are no longer symmetric and the 
resulting solutions are considerably more complicated. Defining /c 12 = k, k 2 ± = ke~ A , and 
A = (ei — e2)/fcsT, we obtain, 



P 2 (t) = l-P 1 (t), (A12) 
^, /u 1 A JU (rArf + Ar t + B) + C) . . . , 

a „ &e~ A 
-Pi + 



and 



2J 2J 



»<**(*)) = ^£ £ ( ; .,-.., fj (,..^.,) ^) 

27o -9f(pi 2 )--(l-2P 1 ) 



where, 

1 / o 2 4 / 3 (a 2 -36) 



ri 



2a + ^ ° U) + 2 2 / 3 t 1/3 , (A15) 



r 2 = r 3 



6 v e 1/3 

./ ^uy-^ ^,,^!, (A16) 
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f = -2a 3 + 9ab - 27c + 3VsV-a 2 b 2 + 46 3 + 4a 3 c - 18abc + 27dc 2 , (A17) 



a = 470 + a, (A18) 

6 = 4J 2 + 4 7o a + |A| 2 , (A19) 

c = 8J 2 7o + a|A| 2 , (A20) 

A = 4 7o + fce" A , (A21) 

B = 4 7o A;e- A + 2 J 2 + |A| 2 , (A22) 

C = 4 7o J 2 + £;e- A |A| 2 , (A23) 

A = icu - 2 7o , (A24) 



and 



a = A;(l + e- A ). (A25) 

pn(i), 5R(/?i2(t)), and ^(pi^i)) are plotted in Fig. El for the parameter set J = 7 o = fc = 
ti = k B T = 1, and e 2 = 0. For unequal site energies the equilibrium values of the classical 
populations do not satisfy the Boltzmann distribution and the coherences are not zero. 

The expectation value of the energy of the dimer is 

(E{t)) = eipn(t) + e 2 p 22 (t) + 2Mp 12 {t), (A26) 

and, as shown by Fig.® although the initial relaxation is faster than in the classical limit the 
equilibrium value exceeds the classical value. Quantum mechanical derealization onto the 
higher energy site increases the average energy, because Pi is larger than its Boltzmann value, 
thus raising the potential energy. This increase in potential energy is partially compensated 
by a reduction in the kinetic energy arising from the last term in Eq. (|A26|) . Fig. El also 
shows the effect of different site energies on (AE(t)). 
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FIG. 5: pn(t) (solid curve), 5ft (/O12 (*)) (dotted curve), and ^(^(i)) (dashed curve) versus time for 
the dimer with J = 70 = k = t\ = k^T = 1, and €2 = 0. Also shown is pu(t) (solid curve with 
crosses) in the classical limit (J = 0). 

APPENDIX B: SOLUTION OF THE EQUATION OF MOTION OF THE DEN- 
SITY MATRIX BY LINEAR ALGEBRA TECHNIQUES 

As described in Section III, for times greater than r the dynamics of the density matrix is 
determined by eigenmodes of A whose eigenvalues satisfy dt\\{\ < r _1 . If only this long-time 
behavior is required it is only necessary to diagonalize A for a small sub-set of the entire 
spectrum using sparse-matrix diagonalization techniques^. Fig. shows Sft(Ai) for a lattice 
of 25 x 25 sites for J = 0.00414 eV. The first 100 eigenvalues satisfy 9ft|Aj| < 9.53 x 10~ 4 eV, 
corresponding to eigenmodes decaying over timescales of t > 4.34 x 10 -12 s. 

Since A is non-hermitian, the inverse of S (where the columns of S are the eigenvectors 
of A) is not equal to its adjoint, and thus its inverse must be calculated explicitly. However, 
if only N out of the total of M eigenvalues and eigenf unctions have been computed, S is a 
M x N matrix, and the matrix equation, 



S - S 



-1 



1 



(Bl) 



is over-determined. Its solution is determined by, 



S 



-1 



V-[diag(l/A 4 )]-Ut 



(B2) 
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FIG. 6: (E(t)) (solid curves) and (AE(t)) (dashed curves) versus time for the dimer with J = jo = 
k = ei = ksT = 1, and e<i = 0. The curves with crosses are the classical limit (J = 0). Notice 
that although the energy relaxation is initially faster in the mixed regime, the equilibrium energy 
exceeds the result of the classical limit. For equal site energies {AE(t)) = \ J\. 

where V and U are determined by the singular-value decomposition of S 21 ' 22 , 



(Note that in practice, as Eq. (fTTj) indicates, to compute {p(t)} only the mth column of S -1 
is required.) 

Fig. IS] compares (r 2 ) 1 / 2 and (E) for 25 x 25 sites obtained by the sparse matrix diagonal- 
ization method (using the 100 highest eigenvalues shown in Fig. |7j) and the Crank- Nicolson 
method. As expected from the eigenvalue spectrum, the results agree for times t > 10~ n s. 



S = U • [diag(Ai)] ■ V\ 



(B3) 
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FIG. 7: 5ft(Aj) versus i for the 100 eigenvalues satisfying 5ft|Aj| < 5R|Aioo|- ^(\) is off-set from zero 
by the non-zero radiative transition rate. 
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t is) t (s) 



FIG. 8: Comparison of (a) (r 2 ) 1 / 2 and (b) (E) obtained by the sparse matrix diagonalization 
method (solid curves) and the Crank-Nicolson method (dotted curves) for a 25 x 25 lattice with 
J = 0.00414 eV. 
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